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Abstract 

We report on a considerable improvement in the technique of measuring multipar- 
ticle correlations via integrals over correlation functions. A modification of measures 
used in the characterization of chaotic dynamical sytems permits fast and flexible 
calculation of factorial moments and cumulants as well as their differential versions. 
Higher order correlation integral measurements even of large multiplicity events such 
as encountered in heavy ion collisons are now feasible. The change from "ordinary" to 
"factorial" powers may have important consequences in other fields such as the study 
of galaxy correlations and Bose-Einstein interferometry. 
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1 Introduction 



Wherever statistical analyses are done, whether in physics, biology or psychology, the mea- 
surement of the correlation function is a basic element of understanding. While each dis- 
cipline has its own set of questions for which it seeks answers, the underlying statistical 
mechanisms are very similar: given a set of variables, one first finds the distribution of 
how often each of these variables takes on a certain value, and much that is useful can be 
learned from these one- variable distributions. Following immediately is the question how two 
variables behave simultaneously, whether they are independent or in some way correlated: 
the two-variable correlation function provides the answers. Higher orders provide additional 
information, but with escalating difficulty of measurement and diminishing returns. 

Conversely, a knowledge of correlation functions to all orders provides complete informa- 
tion on any statistical system. 

Of special interest to us here are point distributions. Typical examples of point distribu- 
tions include galaxies in the sky, cows in a field and particles in phase space (the exact size of 
the object under study is irrelevant as long as it is small compared to the embedding space). 
The aim of this paper is to develop and extend methods of measuring correlation functions 
of point distributions. While we shall be considering particle correlations in high energy 
collisions, the formalism developed here should be suitable, with appropriate modifications, 
for problems in a number of other situations. 

Beyond traditional methods, recent advances in the theory of fractals and scaling have 
spawned a new way of approaching correlations: scaling behavior manifests itself in power- 
law behavior of the correlation function, which in particle physics can be measured experi- 
mentally through the factorial moments revived by Bialas and Peschanski [|l|. On the other 
hand, a measure termed the correlation integral [0, |3], [|] has been in use in the characteri- 
zation of strange attractors and other contexts for some time. As an improvement on the 
factorial moments and following the suggestion of Dremin ||, we previously introduced into 
multiparticle physics two forms of the correlation integral which we termed the "Snake" 
and the symmetrized "GHP" integrals 0. In this paper, we advocate a slightly different 
form, which for obvious reasons we name "Star" integral. Yet another form useful for pion 
interferometry takes the invariant mass of the g-tuple as a measure of its size [0, . 

To illustrate the concept of the correlation integral and the difference between the various 
forms, we consider the phase space plot (e.g. in rapidity and azimuthal angle) of the pions 
produced in a particular collision and ask the question how clusters consisting of q particles 
may be characterized, i.e. what "size" should be assigned to each cluster. Taking q = 5, 
for example, we show in Figure 1 one particular selection of five particles from this event 
and assign the following sizes to it: in (a), the five particles are joined into a "snake", and 
the 5-tuple is assigned a size e corresponding to the longest of the four joining lines. In 
(b), the same 5-tuple is assigned a size given by the maximum of all ten pairwise distances; 
this defines the "GHP" correlation integral. The "Star" integral in Figure 1(c), finally, is 
assigned a size corresponding to the longest of the four lines to the chosen center particle. A 
line represents a particular interparticle distance being tested and found to be smaller than 
e; those particles not connected by a line may or may not be within a distance e of each 
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other. 

Every correlation integral of order q thus assigns a size to every possible g-tuple of 
particles in an event. The way this assignment is done distinguishes the different versions 
of correlation integrals. Once such an assignment is made, they all count the number of 
g-tuples that are smaller than a given size. For a large data sample this corresponds to 
an integration over the q-th order correlation function, as shown explicitly in Eqs. (|4])-(|7|) 
below. 

We shall show that the Star shares all the advantages of the Snake and GHP forms but 
is more amenable to intuitive understanding and is computationally more efficient by orders 
of magnitude. Furthermore, our formalism leads naturally to a conceptual cleanup of the 
heuristic fractal measure used in the study of strange attractors, in galaxy distributions and 
in multiparticle correlations |3|, |], || ||. Being derived directly and rigorously from the 
underlying correlation function, it necessitates the use of factorial powers rather than the 
ordinary powers used previously. The difference becomes significant when small particle 
numbers are involved, a situation which occurs inevitably when the integration domain 
becomes sufficiently small. 

Besides this cleanup, we provide here a technique to measure integrals of cumulant cor- 
relation functions which are the genuine higher order correlations and thus very useful ob- 
servables. Previously, the only way to extract information on higher order cumulants had 
been via combinations of factorial moments [|l(| and in third order for some very special 



configurations in rapidity space [jTl j. By clarifying and extending the technique of event 
mixing, we show how cumulant functions can be integrated directly and thus share all the 
advantages of correlation integrals over bin-averaged factorial moments. 

The paper is organized as follows. In Section |2|, the Star correlation integral is introduced, 
while in Section |3| we lay a solid theoretical foundation for the practice of event mixing. In 
Section |], we remind the reader of the basic relations defining true correlations and show how 
these may be measured directly in the correlation integral scheme. Differential versions and 
their considerable advantages are discussed in Section [5| We conclude with some remarks 
and recommendations. 



2 The correlation integral 

2.1 Basic concepts 

The starting point for all correlation analyses is the g-th order correlation function or "in- 
clusive density" p q (xi, x 2 , . . . , x q ). The choice of the q variables Xk living in d dimensions 
is determined by our particular problem; in high energy physics, we can have, for example, 
vector three-momenta x = (p x ,P y ,Pz), some combination of boost-invariant variables such 
as rapidity, azimuthal angle and transverse momentum (y, 4>,pj_) or just one of these alonef]. 

5 In addition to a suitable choice of variables, some "pre-processing" of the data may be required to 
eliminate unwanted effects, for example the transformation to jet coordinates in e + e~ collisions fll^j , and, if 
desired, the creation of data subsamples of fixed total multiplicity. 
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The set of correlation functions p q , q = 1, 2, . . . can be denned for a fixed total number of 
particles and/or for specific particle types such as positively and negatively charged pions; 
for the purposes of this paper, we shall mostly consider only one type of particle within an 
inclusive distribution where the total number of particles per event is not held constant. For 
this case, p q is defined operationally as 

1 d q a mcl M , 



Pg(x U . . . , Xg) 



<7i dx\ dx2 . . . dXr 



where 07 is the total interaction cross section and <7i nc i the inclusive cross section. Integrating 
in d dimensions over the total window Ax = (Ax) d , we get 



/ 

J AX 



dx 1 ...dx qPq (x 1 ,...,x q ) = (N [q] ) , (2) 



where N is the total number of particles in Ax and (^N^ = (N(N—1) . . . (N—q+1)) = 
£ q (Ax) is the unnormalized q-th order factorial moment over the same region. 

Integrating p q over various domains of integration, one can obtain any number of possible 
moments. For example, the vertical normalized factorial moment revived by Bialas and 
Peschanski Q, 

pv _ _j_ ( n m 1 ,...,m d ) 1 Jn(m) life dx k p q (g], , ■ ■ ■ , X q ) 

q ~ M* m,hn d =i (n mi ,..., md ) q ~ mi ^ d =l /n(m) Uk d^k Pi (*i) ■ ■ ■ Pi (x q ) ' 

integrates p q in a cartesian lattice of M d cubes, each of size Q(m) = Q(mi, . . . , TOd) = {Sx) d , 
normalizing each bin separately. For the purposes of searching for power-law behavior of the 
correlation function, F q is measured as a function of M, with the bin edge length decreas- 
ing correspondingly, 5x = Ax/M. The "horizontal" factorial moments [fl[, the differential 
versions of Section || and indeed the traditional way of presenting correlation functions in 
constant bin sizes are all integrations over different domains of the same correlation function. 

In the same way, the three correlation integrals are simply integrals over specific domains. 
The g-th order Snake correlation integral is defined in terms of p q and the "cluster size" e as 

F Snake/ a _ / Pg( x l, ; ; ; ; x g) Q12Q23 ■ ■ ■ @ g -l,g dx l . . . dx q 
9 J p^CEx) . . . P l{Xq) 613623 • • • ©9-1,, dx t . . . dx q ' 

where 6^ = Q(e — \xi — Xj\). Similarly, the GHP moment is defined with all interparticle 
distances restricted^, 

^GHP / \ _ / Pq( x ii • • • i x q) ri)i<j Qjj dx\ . . . dx q 
q ~ ! pi{xx) .. .pi(x q ) UiKjQijdXi .. .dx q ' 



where 



II % = 0i2e 13 ...e lg e23...e g _ lig (6) 

i<j=l 



6 In Ref. p| , we erroneously denied that the GHP integral could be written down as an analytical integral 
of the correlation function. 
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restricts all possible pairs of coordinates. The Star integral is defined as 



F Star/ x = J PgOgjj LLLL gg) Ql2Ql3 • • • gig ^1 • • • dx q 

~~ fpi(xi) ...pi{x q )® 1% ®xz...® lq dxi...dxq 



involving all coordinates paired with x%. The "topology" of the different correlation integrals 
shown in Figure 1 is thus visible already in the selection of theta functions. 



2.2 The Star integral 

While the definition ([]]) of p q is exact and unambiguous, we shall need for the derivation 
of the correlation integral an alternative but equivalent formulation written down by e.g. 
Klimontovich Ul3 |. Let the N particles of a specific event be situated at the points X{ in 
phase space, i — 1, . . . , N. Then the "event correlation function" p q is defined as adding 1 at 
every point (xi, . . . , x q ) if there is simultaneously at each Xk a particle X, independently of 
the positions of the other N — q particles. This is done for all = N\/(N — q)\ selections 
of g-tuples out of the total N particles. For a specific event a this defines a function 

N 

p«(x l7 ...,x g )= £ 8( Xl - XI) S(x 2 -XI)--- 8(x q - XI 



=1 



l q ) 5 



where S(x) is the product of d one-dimensional delta functions. This function, when averaged 
over all events, yields the g-particle distribution function of Eq. ([[]), 

_ iVev 

Pq (Xt, ...,X q ) = (pgiXi, ...,X q )) = Nj J2 Pq(Xl, ...,X q ), (9) 

a=l 

where A^ev is the number of events in the experimental sample. For finite resolution, Eq.(^[) 
corresponds to building up a gd-dimensional histogram^]. 

Inserting Eqs. into the numerator of Eq. ([/]), we find immediately the (unnormal- 

ized) Star integral factorial moment to be 

ef ar (e) = ( E 8(c - X lll2 )Q(e - X hi3 ) . . . 6(e - X iliq )\ , (10) 

where X^ = \Xi 1 —Xi k \. Pulling out the first sum, we can factorize the remaining sums, 

ef "(e) = (E [ E ®(z - XiiJ) V (ii) 



7 Since both definitions and (^|) are implemented for a finite number of events, they are, strictly 
speaking, estimators of the true correlations. 
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where the factorial power [q—1] in the exponent came about because the sum indices are 
restricted to i a ^ ip ^ i\ for all a ^ [3. The quantity inside the round brackets is so 
important that we give it a special name, the sphere count 

N 

n(X h ,e) = J2 @( e - X nii) , «2^i, (12) 

i2 = l 

which represents the number of particles within a sphere of radius e centered on the particle 
X ix , excluding the center particle itself. (For the given e and center particle shown in Figure 
2(a), we would have n{X^ e) = 9.) In a similar form, it is used extensively in the description 
of galaxy distributions. Introducing the shorthand notation n{X a iv e) = a the unnormalized 
factorial moment can be written compactly as (we henceforth drop the "Star" superscript) 

= a^EE^u,^ 11 = (e« m ) • (is) 



An alternative derivation of £ 9 (e) proceeds by coordinate transformation || |14 |. One 
first defines a "particle-centered" correlation function around the particle at X^, fixing to 
it the coordinate X\ by a delta function, 

N 

p q (X h ;x 1 ,x 2 ,...,x q ) = 8(x 1 -X il ) E 5(x 2 -X i2 )5(x 3 -X h ) ■ ■ ■ 8{x q -X iq ) , (14) 

where the sum indices are all restricted additionally by i a ^ i\, and then transforms to 
relative coordinates. These are the distinctive hallmark of correlation integrals: for the 
Snake integral, we used the coordinate transformation R = J2l=i x k/q. an d = Xk+i — Xk 
0. For the Star integral, on the other hand, all coordinates are defined relative to x±, 

r k = x k -xt, k = 2,...,q. (15) 
Inserting these into the delta functions of Eq. (|T4"D , we find 

p q (X H ;R,r 2 ,...,r q ) = 5[R-X tl ] E 5[r 2 - (X h - X h )\ x . . . 

...x5[r q -(X, lq -X ll )]- (16) 

Since at this point we are only interested in correlations as a function of relative distances, 
we integrate out R over the entire interval AR = Ax. The relative coordinates r k we want 
to restrict to a maximum length r k = \rk\ < e. For higher dimensions d > 1, we must 
first integrate out the angular parts dQ k of dr k . Since we shall eventually normalize our 
correlations using exactly the same domain of integration, however, the constants resulting 
from the angular integrations will cancel and we henceforth ignore them. The remaining 
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integral over the lengths r k is given by J^Ilfe dr^r^ . At the same time, on integrating out 
the angular coordinates, the remaining delta functions of Eq. (|l~6| ) become 



J dVt k 5\r k 



= sln :*r ] . (17) 

Wit 



and the factors 1 will on integration cancel exactly with Xf xi }. To express this entire 
process of simplification, we shall write in shorthand 

l[dr k , (18) 

k 

which is just an integral of the lengths r k . For d = 1, J e Yik dr k is shorthand for e fife dr k . 

Integrating in this way over all relative coordinates, we again get a factorial product of 
single sums: 



/ dR f T[dr k p q (X il ;R,r 2 ,...,r q 
J ax Jo £r 2 



£ e(e - X ili2 )Q(e - X ili3 ) . . . 0(e - X iliq ) 

\ [9-1] 

J2&(e-X Hl2 )) , (19) 



which, on summing over i\ and averaging over all events, again yields Eq. fliTI). 

As derived above, the £ ? (e) of Eqs. (0) and (|TT|)-([T3|) is the unnormalized factorial moment 
over the domain shown in Fig. 2(a), 

fl s (e) = {R, r k \R e Ax, r k e [0, e], k = 2, . . . , q} ; (20) 

this "Floating Sphere" form is best whenever the nature of our variables x in d dimensions 
is such that a length can be sensibly defined, the most obvious being the euclidean distance 



r k = V4* + r L + "-> ( 21 ) 

but often this is not so good. The variables {y,<j),p±), for example, have very different 
behavior, and it may be better to treat each one separately. For such cases, one may use the 
"Floating Box" form [I5|, where each of the d components is treated as a one-dimensional 
correlation integral and the total domain as shown in Figure 2(b) is 

fi fl (e) = {R f , r kJ \R f E Ax, r kJ e [-e, e] Vk = 2, . . . , q; f = y, 0, . . .} . (22) 

This corresponds to inserting d(q— 1) theta functions into Eq. (0), one for each component. 
The results of this section have the following important consequences: 
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1. As shown in Eqs. (|T9|)-(|TT|), the q— 1 sums factorize nicely into a single sum of sphere 
counts n(Xi 1 ,e). This means that £ g (e) can be calculated in an algorithm of order 
N 2 , independently of the order q. As emphasized previously ||, this represents a 
tremendous savings in CPU time over other correlation measures, including the Snake 
and GHP integrals advocated by us previously ||, which run under N q and N q /q\ 
algorithms respectively. 

2. This means that the correlation integral can now be used also for correlation analysis in 
heavy ion collisions, something hitherto impossible due to the large event multiplicities. 
The big improvements in statistics over the conventional vertical factorial moments will 
allow for much more accurate measurements. 

3. Unlike the Star integral found so far in the literature || ^, 16| , we get a factorial 



power of the sphere count n iq ^ 1 ^ rather than the ordinary power n q ~ 1 . This result we 
obtained rigorously from first principles merely by restricting the sum indices to be 
unequal because the same particle may not be counted more than once. Even for the 
large multiplicities encountered in astronomy, the change is not inconsequential, as the 
difference between n iq ~ 1 ^ and n q ~ x is important when e becomes sufficiently small. 

4. We obtained this factorial power without drawing on distinctions between "dynamical" 
and "statistical" fluctuations 

5. To first order, we have ignored the variation of £ 9 (e) with the center coordinate R; 
this is equivalent to assuming that the physics is the same for all parts of the defined 
window. When statistics permit, it may be very useful to measure £ g also as a function 
of R. For example, one expects the correlation structure at small transverse momentum 
to be very different from that at large p±, so that a separate measurement of £ 9 for 
small and large R = p±\ may be most revealing. 

6. The advantages of correlation integrals over the traditional factorial moments arise 
because the former use interparticle distances directly while the latter rely on fixed 
bins and grids rather than the particle positions themselves. 

7. We must emphasize that the Star integral as derived here is different from the Snake 
and GHP versions we used in earlier papers. All three are correlation integrals, but 
they differ in the topology of the interparticle distances measured. For q — 2, all three 
are the same, while the Snake and Star versions are the same even at the level q = 3. 
Only in fourth order do the differences between the latter two appear. While the Snake 
and GHP versions are not wrong and an improvement over previous work, the present 
Star integral represents a big step forward. 
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3 Normalization by event mixing 



We now proceed to consider the denominator of Eq. (^). Since the topic of normalization 
is complex and full of pitfalls, we do not address the full range of issues here and defer 
such discussion to future work. Instead, we concentrate on the normalization scheme we 
consider most suitable for the measurement of correlations in high energy physics, a version 
resembling the so-called vertical normalization used in Eq. ([3D. 

The normalization of the correlation integral is done by means of event mixing, a seem- 
ingly heuristic technique commonly used in Bose-Einstein correlations ||12||. It is well founded, 



however, both for our purposes here and in the Bose-Einstein context ||. We recall that if 
all q coordinates x k are statistically independent of each other, the correlation function fac- 
torizes: p q = p\. We hence normalize the numerator (|TT| ) of the Star integral by integrating 
llfe=i Pi( x k = R+fk) over the same domain, 



dR Pl {R) Y[dr kPl (R+r 2 ) ... Pl (R+r q 

AX Jo ^ = A 2 



AX 



dR Pl (R) 



dr 2 pi{R+r 2 ) 



q-l 



(23) 



(here and below it is understood that one transforms from x- to r-coordinates before integra- 
tion). Inserting p\(x k ) = N^ 1 J2e k J2i k $( x k — ^-tt) f° r eacri factor, we find after integration 



e orm (e) 



ac 1 EE 



e2,...,e 9 i 2 ,...,i„ 



(24) 



where now Xf}? h = 



I \^ e l 



Y e k I 



measures the distance between two particles taken from 
different events e\ and e k - This much resembles the numerator expressions of Eq. flTTD, and 
indeed the sums also factorize here, so that 




E(^W,e)) 



9-1 



(25) 



or, in shorthand, 



™ (e) 



£< 6 > 4 



(26) 



Comparing the numerator ( 
latter is an ordinary power q 
from the fact that in Eq. (p 
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to the denominator fl2q) , we note that the exponent of the 
1 instead of the factorial power [q— 1] of the former; this follows 
there are no restrictions on the sum indices. 
Figure 3 shows how the denominator sphere count n(,(X", e) = b is found by inserting the 
particle X.\ of event a into another event b of the sample and doing the count n b (X^, e) = 
J2j ©(e ~ around it. This is done for many events b to obtain the inner event average 



of Eq. fl25|). Note that one has to distinguish carefully between a- and 6-event averages: 
the computation of (b) = (n b (X",e)) involves an average over different events b while the 
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position of the sphere center is kept fixed at X". Expressions of the form (J2i(b)) then denote 
sums of contributions when the center of the sphere "floats" over the i = 1, . . . , N particle 
positions of event a and finally over all events a — 1, ... , N ev . 

While the derivation of Eq. (|25|) from (|23D is exact, this expression for the normalization 
is correct only for N cv — > oo: it contains a hidden bias due to correlations induced by use of 
the same events in every factor pi in (^). For the case where the the experimental sample 
size is not infinite, statistical theory provides estimators, which from the limited-size sample 
estimate quantities for the "true" (i.e. infinitely large) sample. Applying the theory of 
estimators to our problem, we find that the product of distributions p\ in Eq. (|23| ) must be 
modified precisely in such a way that the event indices all mutually unequal. 

This is in agreement with the heuristic procedure of creating "fake events" where each track 
is taken from a different (real) event. 

The corrections to obtain such an unbiased form of the normalization can be written as 
a series in powers of l/N eY , with the leading term given by For the relatively small 

number of events and great sensitivity to bias found in heavy ion samples, the correction can 
be quite important, while the situation for hadronic data is less acute. We defer the details 
of this very technical discussion and the exact expressions to a future publication |17|] . 



Subject to the above corrections, the normalized q-th order moment in its Star integral 
form is 

F(£) = JifcL J£' a "-") (27) 

As in the case of the traditional moments (|3|), one would measure F q (e) as a function of 
decreasing e; a straight line in a plot of \nF q vs. lne would, as before, be interpreted as 
scaling behavior of p q . We caution, however, that there are important issues which must be 
addressed before such conclusions can be drawn, among them normalization effects. 
The following points should be noted. 

1. Theoretical models are easily compared to data obtained with the Star integral: when 
they predict the single-particle distribution pi and correlation functions p q , their cor- 
responding Star integral moment is just given by the analytic integral expression (0). 
Monte Carlo simulations should, of course, take into account the proper removal of 
bias induced by finite sample size. 

2. The measurement of the Star integral is very economical. Just as the numerator £ g can 
be measured in an algorithm of order N 2 (the number of particles per event), so the 
denominator requires only an algorithm of the order of the square of the sample size, 
N 2 V , a large savings over the previous GHP and Snake integral algorithms 0. (This 
savings is not destroyed by the abovementioned corrections to obtain an unbiased 
normalization.) When the order of the events in the sample is random, this can be 
reduced even further by taking for the "inner" event average (6) only a fraction of the 
full sample, e.g. b G {a— I, a— 2, . . . , a— A}. The smaller A, the faster the calculation 
but the larger the statistical error. Since a small value of A introduces a considerable 
bias, which disappears as it is increased, great care must be taken that A is of a size 
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where the normalization becomes independent of its exact value. An optimal value for 
A can be found for a given sample and length e by experimentation. Whenever doubt 
arises, the full A = N ev —1 sum should be taken as this is exact. 

For smaller values of A, the error in the denominator will be nonnegligible and should 
be combined with the numerator error, including covariances. 

3. For e = Ax, all theta functions become unity by default and so F q (Ax) = (N [q] )/(N) q , 
which is unity only when the event multiplicities in Ax are Poisson-distributed. Fur- 
thermore, when the distribution is purely random for a given e, then F q (e) becomes 
unity; see the discussion around Eq. (|50|) below. 

4. Annoying boundary effects due to the fmiteness of Ax are largely cancelled out, because 
the vertical normalization used means that sphere counts for centers Xi close to the 
boundary are reduced for both numerator and denominator. In the business of galaxy 
distributions, this remains a much-discussed topic 0; since unfortunately only a single 
event exists in this case, the horizontal normalization, which is vulnerable to such 
boundary effects, has to be used there. 

5. Nevertheless, there may be instances where the horizontal normalization may be pre- 
ferred. The definition of fractal dimensions, for example, is often couched in terms of 
the horizontally-normalized correlation integral. 

6. The vertical normalization has the additional advantage that the effects of the single- 
particle distribution pi(x) are cancelled out to some extent. While in many cases 
pi is constant ("stationary") or varying only weakly, it may in high energy collisions 
vary quite strongly; especially the transverse momentum distribution p\(p±) is strongly 
peaked and then falls off exponentially. Since sphere counts in both numerator and de- 
nominator vary about equally as a function of pi, the trivial effects on F q of a strongly 
varying single particle distribution are compensated. However, this compensation is 
only partial: when e is small, each small domain is approximately flat and the cancel- 
lation is fairly reliable. When e is large, though, the variation of a nonstationary p x 
within the domain is averaged out before the ratio is taken instead of the better reverse 
order. Ideally, one would divide up such large domains into many small ones and sum 
up the contributions only after normalization, but this is usually made impossible by 
bad statistics. Also, even the so-called vertical factorial moments (0) used traditionally 
suffer for large bins from the same effect, so that this problem can be ignored at the 
present level of sophistication. 

7. Of greater concern is the possibility that the integration of the center-coordinate R 
over the entire domain Ax in Eqs. (Tl9| ) and fl23| ) suffers from the same problem for 
any value of e. For pi's that are weakly varying this does not matter very much, but 
when they change drastically (such as Pi(p_l)) a more vertical form should be used. At 
some cost in CPU time, this can be implemented straightforwardly. Instead of letting 
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R range over the entire space, we introduce discrete binning; for the Floating Box 



By splitting the R integration in this way, the correlation integral can be made largely 
independent of the shape of p\. As one is usually interested in F q as a function of e 
only, the number of subdivisions L and bin sizes SR can be kept fixed throughout a 
computation; the only requirement is that the number of subdivisions should be large 
enough to ensure that the final correlation integral depends on L only weakly. 

For systems where p\ varies strongly only with one of its variables while varying weakly 
with the others, it may be possible to define a "hybrid" correlation integral with mixed 
normalization, implementing an .R-sum as in Eq. (|28|) for this component only and not 
for the weakly-varying components. This would be useful to boost the number of 
counts. 

As discussed in Section |272|, splitting up R into different regions may in itself be useful 
in isolating different physical effects. No sums over i would be taken for such cases. 

8. There remains the question how the traditional factorial moments (§) and Star corre- 
lation integrals ( $F7\) are to be compared. The simplest answer is that they should not 
be compared at all, since by their different definitions they cannot be expected to yield 
exactly the same results even for the same data sample. Under ideal circumstances 
(N ev — KX3, weakly varying p 1; etc.), the two should yield a similar slope asymptoti- 
cally. It may be helpful to compare the two at equal e, but deviations should not be 
interpreted as revealing anything fundamental. If one finds large differences between 
traditional factorial moments (§) and correlation integrals, then they are mostly due 
to the different normalization procedures used. 

One can, however, check for consistency between the BP factorial moments and the 
Star integrals. Both should be the same for e = Aa;. Also, for integer M, one can do 
the Star integral separately for every bin, with e set to Ax/M, and then add up the 
contributions. This sum over Star integrals should then be identical to the vertical 
factorial moment F^(M). 




(28) 
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4 Cumulants 



4.1 Definition and use 

The measurement of factorial moments and/or correlation integrals may be useful in itself in 
searching for a power law in the correlation function. Moments do not, however, reveal the 
true correlations, because the correlation function contains uncorrelated parts which have 
to be subtracted. This becomes clear on considering the effect that statistical independence 
has on the correlation function. Statistical independence of two points X\ and x 2 means that 
the correlation function p 2 (xi,x 2 ) factorizes into a product pi{x\) pi(x 2 ) . Similarly, when 
x\ becomes statistically independent of all other points Xk^i, the higher order correlation 
functions factorize accordingly: p q (x 1 , x 2 , ■ ■ ■ , x q ) — > p 1 (xi)p q ^i(x 2 , . . . , x g ). All possible 
combinations of such factorizations have to be subtracted from the original correlation func- 
tion before one can speak of the "true" correlations. 

These reduced quantities, known as cumulants, are basic to statistical analysis of any 
sort |18| . They are constructed precisely in such a way as to become zero whenever any one 



or more of the points becomes statistically independent of the others. (The often-used 
factorization p q — > p\ is only the most drastic form of statistical independence, assuming 
that every point becomes independent of every other.) Cumulants of different distributions 
are also additive under convolution of the distributions |TjJ as well as being invariant under 



change of origin [18 



The first few cumulants C q are, in terms of the correlation functions, 

C 2 (xi,x 2 ) = p 2 (x 1 ,x 2 ) - pi(xi)pi(x 2 ) , (30) 
C 3 (x u x 2 ,x 3 ) = p 3 (xi,x 2 ,x 3 ) - px {xx)p 2 {x 2 , x 3 ) 

- Pi(x 2 )p2(x 3 ,X 1 ) - Pl{x 3 )p 2 {x U X 2 ) 

+ 2pi(aji)pi(aj 2 )pi(aj 3 ) , (31) 
C^x^x^x^xt) = p i (x 1 ,x 2 ,x 3 ,x 4 ) - J2pi( x i)p3(x 2 ,x 3 ,x 4 ) 

(4) 

- ^p 2 (xi,x 2 )p 2 (x 3 ,x 4 ) 

(3) 

+ 2 Pi ( x i)Pi (x 2 )p 2 (x 3 ,x 4 ) 

(6) 

- Gpi{xi)p 1 (x 2 )p 1 (x 3 )p 1 (x 4 ) . (32) 

The bracketed numbers under the sums indicate the number of permutations of the argu- 
ments Xk which have to be included. Further, omitting the arguments, 

C5 = P5 ~ J2 PiP* ~ J2 PiPz 

(5) (10) 

+ 2 J2 P1P1P3 + 2 J2 P1P2P2 

(10) (15) 

-6j2piPiPip2 + 2Apl, (33) 

(10) 
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and so on for higher orders. 

These equations have been utilized to find simple relations between the vertical factorial 
moments F q of Eq.(|3]) and the integrated normalized cumulants [1C]. With m = (mi, . . . , ma) 
as usual, the integrated normalized cumulant is defined as 



K{M) 



1 



E 



In m Ilk dx k C q (xi, . . . , x g ) 



Sa m dx Pi( x ) 



(34) 



which yields relations such as K 2 = F% — 1, K% = F% —3F 2 +2 etc. which can thus be utilized 
directly by experimentalists. These relations hold exactly for q < 3 and approximately for 
q > 3. They are not true for horizontally normalized moments. 



4.2 Correlation integral cumulants 



In contrast to Eq. (p4|) , the correlation integral cumulant is defined as the integral of C q over 
the domains Q$ of Eq. (|20| ) or of Eq. (|22| ) after appropriate transformation to relative 
coordinates, 

K q (e) = - ft?(£ ) . , (35) 

with 

f,(e) = / dR f e l[dr k C q (R,r 2 ,...,r q ) 

J AX JO u 

= J C q (xi, ...,x q ) 612613 . . . Q lq dxi . . . dx q (36) 

the unnormalized factorial cumulant. The latter can be written entirely in terms of the 
sphere counts introduced previously, 

j 

b = ^e(e-X^)=n 6 (X«,e); (37) 

j 

(see Figures 2 and 3 for terms contributing to a and h respectively). To demonstrate this, 
we start with q = 2. Here f 2 = f(p 2 — PiPi), which by Eqs. (|TTJ) QT^) and fl2"3|)-(PBD is seen 
to yield (henceforth we suppress the dependence on e), 

h = foa - (b))\ . (38) 

For q = 3, the first term in the expansion fl3T|) of C3 is just £3 = 0^), while the last term 
P1P1P1 yields ^3° rm = (J2i{b) 2 )- The three "mixed terms" involving both p 2 and p\ must be 
worked out explicitly. 6n the one hand, 



[dR f e dr 2 dr 3 p 1 (x 1 )p 2 (x 2 ,x 3 ) = N~ 2 £ £ £ 6(e - X«f)6(e - X? k b ) , 

J J0 a,b i j^k 

= (z<& [2] >y (39) 
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while on the other hand, if X\ is contained in p2, 

J dR dr 2 dr 3 p 1 (x 2 )p2(x 1 ,x 3 ) = (j2a(b)^ , (40) 
so that, putting all the pieces together, 

h=(j2(* [2] -(b [2] }-2a(b} + 2(by)\. (41) 



The constant recurrence of the outer event average and z-sum suggests that we define an 
"i-particle cumulant" by 

(E^(o)=^. (42) 

in terms of which we find 





= a - (b) , 


(43) 




= a [2] - (b [2] ) - 2a(b) + 2(b) 2 , 


(44) 




= a^-(b^)-3a^(b)-3a(b [2] ) 






+ Q(b)(b [2] ) +Qa{b} 2 -Q{b} 3 , 


(45) 




= a W_ (6 [4] ) _ 4a [3] (6) _ 4a(6 [3] ) 





- 6aM(bM) + 8(b) (6®) + 12a®(b) 2 + 6(b^)(b^) 

+ 24a(b)(b [2] ) - 36{&) 2 (fe [2] > -24a(6) 3 + 24(6) 4 , (46) 

which are then summed over all particles % and averaged over all events to yield i q . 

How these sums can be obtained graphically is illustrated for q = 3 in Figure 4(a) and 
q = 4 in Figure 4(b). The black squares represent individual particles; those enclosed by a 
circle belong to the same event. The center particle at _X"" is connected to the q—1 other 
particles by the lines representing the theta functions. One now draws all possible event 
topologies with q—1 lines connected to one center particle. For p joining lines within the 
X® event, one writes down a factor a' p '; while lines connecting X" to p particles in the 
same (other) event yields a factor (jffl). p lines going to different events 6, c, . . . results in a 
factor (b) p . Putting all such factors together and assigning to each the appropriate sign and 
prefactor from Eqs. fl3Cf) jf.) one obtains the cumulant expansions (|43|) - (|4~6l) . 

Writing higher orders recursively in terms of lower-order cumulants, 



(47) 



Hi) 


= a®-(b®) 


- 2<6)f 3 (0 


1 


f 4 (z) 


= a®-(b®) 


- mh(i) 


- 3(&[ 2 l)f 2 (z) , 


hii) 


= M-<& [4] > 


- Mb)U(i) 


-6(^)f 3 (i) -4<&[ 3 l)f 2 (i 



we are led to the conjecture that 

9-1 



f ? (i) = - (b^) - Y, [ q _ \) (b [ «- p] )f p (i) , (46 



p=2 
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which, if proven for arbitrary q > 5, could open the way for an easy calculation of cumulants 
to arbitrary order without doing the messy algebra involved. 

Just as the normalization ££ orm must be corrected for bias, we must for third and higher 
order cumulants also correct for limited sample size. Here, too, using estimators to correct for 
the effect of the limited-size sample requires that all event sums go over unequal events. This 
again results in corrections of the order l/iV ev , l/iVf v , ... so that, for example, f 3 (z) above 
would acquire the additional correction term 2((6) 2 — (b 2 )) / (N cv — 2). While these corrections 
can be quite important, we defer a discussion of their origins and exact expressions to future 



work [17 



When an event sample has no correlations, the count around a particle in event a would 
on average be the same as when it were inserted into other events, 



a 



[9-1] 



(b [q - l] ) (49) 



which results in f q = 0. This would of course be true on average only. Less obvious but 
also true is that if any one of the q variables of C q becomes independent, the integrated 
cumulants f q become zero also; this can be shown graphically too. A stronger condition of 
randomness, comparable to the Poisson distribution in fixed bins, would be reflected by the 
behavior 

i i 

(E<^ 1] » — (E^r 1 )- (so) 

i i 

so that the moments would go to unity, F q — > 1, for this case. 

It must be stressed that the g-th order cumulant contains no correlations of order lower 
than q. Thus even if f 2 > 0, f 3 can still be zero when there are no true third order correlations; 
Eq. ([TJ) is merely a convenient grouping of the terms. 

We further see that, since only the basic quantities a and b are needed to construct 
cumulants, they can also be calculated very economically with order N 2 algorithms. 

The estimation of errors always provides a headache since mostly one has to deal with the 
intricacies of error propagation and covariances. In the Star integral formulation, however, 
this process is much simplified; for the statistical error on K q , one merely has to calculate 

, HQ J^fhl (51 , 

1- "ev 

and combine this with the corresponding statistical error for the event mixing denominator. 
Errors for the factorial moments F q and differentials of Section are obtained with the same 
ease. 
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4.3 Presenting cumulants 



While there may be many useful ways to plot correlations, depending on what one is looking 
for, we recommend the following format for cumulants. The second order cumulant 



Me) = TT^V ~ 1 (52) 



cannot be smaller than —1. Its theoretical maximum is harder to find; but an estimate can 
be made by using the extreme case where the events are extremely "spiky" but the spikes 
are uniformly distributed in phase space (from event to event). The event-averaged count 
will be approximately n = (N)(e/Ax) d , while the theta functions of the numerator are all 
unity for our narrow spikes; this extreme case thus gives an approximate upper limit 

/N (N(N-l)) fAx\ 2d , N 

-KK 2 (e)<^- -f-l~ . (53) 



(n(n — 1)) 

to which a given sample K 2 can be compared. For higher orders, one may then plot the 
K q directly and/or as the ratios K q /K 2 , which would express q-th order correlations as a 
fraction of second-order correlations. Testing for linking, on the other hand [|l(J , one would 
plot the ratios Kg/K^ 1 . 

Apart from calculating cumulants which are averaged over the entire event sample, it 
may in specific cases be interesting to look at "single-event cumulants" f q for rare events, 
for example if a certain event or group of events exhibits unusual patterns in a phase space 
plot. To see whether such cumulants differ significantly from mere statistical noise in the 
fluctuations, they should be plotted on top of the corresponding event-averaged cumulant 

plus/minus twice (or three times) the error f q ± 2 f$(*)] 2 ) ~~ fg) ^ • Whether such 
single-event cumulants are a good idea will have to be established in practice. We remind 
the reader that the distribution of factorial moments is not necessarily gaussian [2(J and 
that a large-deviation analysis may also be appropriate in such cases |2l] . 



5 Differential versions 

Correlation integrals and their cumulants described so far are defined always in terms of a 
maximum distance e; the ubiquitous theta functions ensure that all interparticle distances 
Xij involved must be smaller than e. The simplicity of this definition allows one to test 
clusters of many particles at once, i.e. probe correlations of order 3 or higher, something not 
possible using the conventional methods of measuring correlation functions. It makes good 
sense, however, to ask not only whether some interparticle distance is smaller than some 
value, < e, but also whether it falls within a certain distance interval [et-i,et]. 

To this purpose, we define a sequence of distances ei, £2, ... up to some maximum distance 
determined by the total domain of integration. This sequence can be either linear, 

e t = te x , (54) 
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or exponential, 

e t = eic* -1 , c> 1, (55) 

the second definition being useful because data presented as function of In e t will be equally 
spaced. These two sequences divide up the whole phase space into adjacent and disjoint 
domains; for d = 1 and q — 2, these domains are shown in Figure 5 as the sequence of strips 
filling the entire domain. We also introduce the indicator function 

I t (X) = 9(e t - X) - 0(e t _ 1 - X) , (56) 

which is unity when e t _i < X < e t and zero otherwise. 

The differential forms are defined as follows (see also Figure 6). Given a center particle 
X° in event a, the number of particles situated a distance X ilik 6 [e t _i,e t ] away from X i± 
is 

A| 2 (z,t)=7i(Xf 1 ,e t )-n(X? 1 ,e t _i)EEo t -a t _i, (57) 

the latter defining the shortened notation we shall be using. We next ask how many clus- 
ters of q—1 particles exist for which the maximum distance to X^ is in this interval, 
max(Xj li2 , . . . , X iliq ) G [e$_i,e t ]. The answer is surprisingly simple: the number of such 
clusters is 

A&(i,t) = n(Xt^t) [q - l] -n(Xt,e t ^)^ ee af^ - afc 1] , (58) 
since through use of 

0(er-X) = £W (59) 
t=i 

and 

0(e - max(X 1 , ...,X q )) = f[ 9(e - X p ) (60) 

P =i 

we can show that 

/ \ [3-1] / \ fe-l] 

A| ? (ii,t) = (EO^-^mJJ - f£9(e i „ 1 -X M; 

= H J 4 (max(X ili2 , . . . ,X iliq )) . (61) 

In Figure 6, one such cluster is shown, where at least one particle besides the central one is 
in the shaded regions for the Floating Sphere domain of Eq. ( |20| ) . 

The normalization proceeds along the same lines as in Section ^|, so that not surprisingly 
we find that the normalized differential moment is 

Af - (t) = oW- - (*-.>«-r (62) 

while the cumulants are similarly 

(E,f,(it)-f,(M-l)) 
A/ ^ (i) = <E,<M-> - <««-><-> ' <63) 
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Both AF q and AK q are thus accessible to measurement with minimal additional effort. 
Usually, they will be plotted as a function of t, i.e. the distance interval within which the 
maximum interparticle distance would fall. AF q can be thought of as the analogue to the 



"factorial correlator" defined in Ref. [jj and AK q to their cumulants p2 



We conclude this section by listing the advantages of using the differential moments and 
cumulants as a measure of cluster size. 

• First of all, they have all the advantages of the correlation integrals in that their statis- 
tics will be higher and the data points calculated more stable than the corresponding 
correlator. Especially for higher dimensions, the large gain in statistics will permit 
measurements that would be impossible otherwise. 

• Differential moments and cumulants are almost immune to the problems of converting 
from biased to unbiased estimators in the normalization and in f q . Because these 
corrections for bias takes the form of an additive series, taking the difference f g (i,t) — 
f q (i,t—l) causes them to cancel to a large degree. 

• When a unique distance can be defined as in Eq. ([H]), the differential count is unam- 
biguous even in higher dimensions, unlike the correlator, where distance definitions are 
ambiguous for multidimensional analyses 



Very importantly, the domains of integration of the normalized differential moments 
and cumulants are disjoint (Fig. 5), meaning that the data points will not be correlated 
amongst themselves, a constant bug in ordinary moments and cumulants. 

A special status must be accorded to AF 2 (t): it behaves as a "roaming distance", 
looking for all particle pairs that are a certain distance apart. This is suggestive of 



interpreting AF 2 as a kind of Fourier transform of the distances [24 



6 Conclusion 

We have developed a general formalism for measuring correlations of point distributions. 
The language used has been that of high energy physics, but we believe that the method 
may be of use in other fields also. The use of a delta function notation has enabled us to 
derive the Star integral from first principles and through its greater clarity pointed to a 
number of important extensions. Most important of these is that the correlation integral 
used in the astronomy literature || [| and suggested for high energy physics |J appears to 
be in need of modification from the form n q ^ 1 to the "factorial power" form fv- q ~ 1 \ To assess 
whether such a modification is possible or practical in galaxy distributions is beyond the 
scope of this paper; for the limited number of particles encountered in high energy physics, 
however, it seems an unavoidable and clearly superior formulation. 

The felicitous definition of relative coordinates leading to the Star integral makes the 
latter very economical to calculate, including all cumulants and differential quantities. Since 
in addition the domains of integration are the largest possible, we believe that the Star 
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integral extracts the maximum available correlation information from the data for the min- 
imum price in CPU time. In this respect, it has proven superior to the traditional factorial 
moments and older correlation measurements, especially when correlations are measured in 
higher-dimensional spaces. As computers continue their evolution to previously unimagin- 
able speeds and capacity, the present method should become routine even for the higher 
orders. 

Three issues have been dealt with only cursorily or not at all: the question of eliminating 
the influence of the total multiplicity distribution, the measurement of correlations between 
different species of particles (e.g. distinguished by their charge), and the problem of properly 
defining fractal dimensions for nonstationary ensembles of events rather than the usual single- 
event time series. We hope to return to the latter in future. 

The present paper has emphasized the how to rather than the what of measurement. 
At first sight, it may seem unnecessary to devote so much effort to the mere process of 
measurement. However, high energy hadronic collisions, and to an even greater extent 
nucleus-nucleus collisions, are presently in a state where very few exact calculations on 
correlations can be done and most theoretical work relies on assumptions which often are 
hard to support or believe. In this context, we believe that it is of cardinal importance that 
there should be a clear and unambiguous method for analyzing correlations, and that, if 
possible, a standard should be established by which different experiments can be compared. 

The confusion underlying the dynamics of correlations is reflected in the cursory way in 
which we have treated the choice of variables. While there are some theoretical preferences 



25| , we are fairly ignorant of the dynamics of the soft component and hence the best choice 
of variables in this case. The original proposal regarding self-similarity in particle production 
|J did not have theoretical grounding in currently acceptable physical theories but was based 
on a toy model to illustrate the idea, while exact calculations of correlations e.g. in a QCD 
framework [p6] can be applied only to high-p^ processes. 

The occurrence of different forms of correlation integrals might cause unpleasant con- 
fusion among experimentalists who would prefer a unique recipe to extract information of 
higher order correlation functions. Unfortunately, there is a priori no best choice. The dif- 
ferent forms merely reflect the freedom of choice of the particular shape of the integration 
domain; but all commonly probe the correlation functions by decreasing the size of the 
integration domain. While the numerical values of the various integrals may differ, the func- 
tional dependence is supposed to be similar (this was shown numerically for the Snake and 
GHP integrals Moreover, by suitable normalization most of the numerical differences 
between the various forms can be divided out, so that the choice of a particular form can be 
guided by practical efficiency arguments. 

On very general grounds, the choice of relative coordinates seems a wise one, in whatever 
variables one may prefer. This is true especially for cases where there is some degree of sta- 
tionarity in the distribution (i.e. invariance under translation), which is generally assumed 
to be true for galaxy distributions and perhaps for pion distributions at higher collision en- 
ergies. The ubiquitous use of such stationarity assumptions testifies to their popularity. Not 
least, measurements in Bose-Einstein interferometry rely on relative coordinates, whether in 
their three-vector form q = p 1 — p 2 or as a function of Q 2 = — (pi — P2) 2 '■ 
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In this context, we have aimed to provide a framework that is adaptable to any future 
choice of variables and dynamical theories. This will hopefully allow for clean measurements 
to guide theoretical thinking, while remaining flexible in its implementation. 
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Figure captions 



FIG. 1. Different version of the correlation integral: (a) Snake, (b) GHP, (c) Star. For a 
given set of q particles in phase space (here q = 5) taken from the N particles of a particular 
event, pairwise distances are tested according to the topology of joining lines shown. The 
longest of the lines characterizes the size of the given g-tuple in every prescription. The Star 
count is much more efficient than the other prescriptions, see Section ^|. 



FIG. 2. The conceptual advantage of the Star integral over other versions stems from the 
fact that counting g-stars can be reduced to computing factorials of sphere counts n(Xi, e). 
(a) "Floating Sphere": n(Xi, e) = number of neighboring particles within a sphere of radius 
e centered at particle i with coordinates Xi. The center particle itself is not included in the 
sphere count, (b) When the coordinates x have very different physical properties in their 
components (such as y, <p and p±), the "Floating Box" may be a better choice as it treats 
the distances along different coordinate axes independently. 



FIG. 3. The basic building blocks for computing the normalization of correlation integrals 
(Eqs. (j23D-([26|)) as well as cumulants Eqs. fl37|) jf.) is the sphere count hb(Xl,e). While 
similar to the count of Fig. 2(a), n^X^e) performs a sphere count around particle X" 
taken from event a (shown as a dot) by placing it in the event b and counting the 6-particles 
(shown as crosses) within the sphere. In the example shown, h^X 1 *, e) = 6. Averages over 
many events b are taken while X a { is kept fixed. 



FIG. 4. Schematic representation of event mixing terms entering the (a) third (b) fourth 



order cumulant (Eqs. (f4^-([45|)). For a given center particle % in event a, the other particles 
in the g-tuple can be either within the same event a or in different events b,c, . . . (see text). 
p particles in event a lead to a factor a' lP ~ 1 ^; p particles in event b give a factor (b^) and a 
particle in p different events b, c, . . . , gives a factor {b) p . With the appropriate combinatorial 
prefactors, the sum of these terms yields the third and fourth order cumulant integrals. 



FIG. 5. The exponential sequence of distances e t (t = 1,2,.. .) of Eq. (|55D , used to define 
the differential forms of correlation integrals of Section ^. The shaded regions represent the 
integration area of the differential integral AF 2 (t) over the two-particle density p(xi,x 2 ). 
Note that the coordinates x in this figure are one-dimensional (d = 1) and the labels 1 and 
2 refer to two different particles within the interval Ax. 



FIG. 6. An example of the differential sphere count within the shaded area in two dimensions 
(d = 2) gives A£ 2 (M) — a t — a>t-i = 9 — 6 = 3 (Eq. (|5T|)). For higher orders, the number of 
all g-stars with size within the interval [e t -i, e t ], i.e. with at least one of the q— 1 neighboring 
particles within the shaded area, is given by A£ q (i,t) = af~^ — afT^ of Eq. (|58]). In the 
present figure, for example, A£ 3 (i,t) = 42. 
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